function [LE1, LE2] = LEs(func, a, b, c, d, x0, y0, N, discard)
% compute_LEs 计算给定模型 func 下的两个李雅普诺夫指数
%
% 模型为 MNFM1：
%   x(k+1) = a*x(k)*(1 - x(k)) + b*sin(c*cos(y(k))*x(k)) + d
%   y(k+1) = y(k) + x(k)
%
% 输入：
%   func   - 模型函数句柄，如 f1 (应返回 [x_next, y_next])
%   a,b,c,d- 模型参数
%   x0, y0 - 初始条件
%   N      - 总迭代步数
%   discard- 丢弃的暂态步数
%
% 输出：
%   LE1, LE2 - 两个李雅普诺夫指数

    % 初始化状态
    x = x0;
    y = y0;
    
    % 初始化两个扰动向量
    e1 = [1; 0];
    e2 = [0; 1];
    
    % 用于累积扰动向量对数增长率
    sumLog1 = 0;
    sumLog2 = 0;
    
    for k = 1 : N
        % 计算雅可比矩阵：使用通用 jacobian 函数，传入模型函数 func
        J = jacobian(func, x, y, a, b, c, d);
        
        % 更新扰动向量
        e1 = J * e1;
        e2 = J * e2;
        
        % 对 e1, e2 进行 Gram-Schmidt 正交化
        norm_e1 = norm(e1);
        e1 = e1 / norm_e1;
        e2 = e2 - (dot(e2, e1))*e1;
        norm_e2 = norm(e2);
        e2 = e2 / norm_e2;
        
        % 从丢弃暂态之后累计对数增长率
        if k > discard
            sumLog1 = sumLog1 + log(norm_e1);
            sumLog2 = sumLog2 + log(norm_e2);
        end
        
        % 更新状态调用模型函数 func
        [x, y] = func(x, y, a, b, c, d);
    end
    
    M = N - discard;
    LE1 = sumLog1 / M;
    LE2 = sumLog2 / M;
end
